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The simulation of granular particles in a quasi-two-dimensional container under vertical 
vibration as an experimental accessible model for granular gases is performed. The velocity 
distribution function obeys an exponential-like function during the vibration and deviates 
from the exponential function in free-cooling states. It is confirmed that this exponential- 
like distribution function is produced by Coulomb friction force. A Langevin equation with 
Coulomb friction is proposed to describe the motion of such a system. 
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Granular physics has been a challenging field in statistical physics since the rediscovery of 
significant nature in granular materials in the late '80s or early '90s. 1 An assembly of grains 
has strong fluctuations in the configuration and the motion such that mean-field theories 
cannot be used in most situations. A typical example of the strong fluctuation appears in 
the force distribution for a static granular assembly piled by gravity. 2-4 The force propagates 
along force chains and the distribution function of the magnitude of the force on the bottom 
of a container does not obey a Gaussian function but an exponential function. 2-4 

Such a strong fluctuation originaling from non-Gaussian properties should be relevant 
even in the dynamics of granular assemblies. However, there have not been many system- 
atic research studies that focus on the statistical distribution functions in the steady state, 
because real and numerical experiments report no unified results, i.e., velocity distribution 
functions (VDFs) obey Gaussian-like functions with exponential tails, 5_12 from the Gaussian 
to the exponential depending on density, 13 ' 14 stretched exponential 15-17 and even power-law 
functions. 18-21 

As in standard statistical mechanics, an assembly of grains in a gas phase is an idealistic 
system for studing statistical weight. Approximate granular gases can be obtained by rapid 
granular flows on inclined slopes, 22,23 gas-solid mixtures 24, 25 and grains under the external vi- 
bration. 11-18 However, these systems cannot be regarded as idealistic granular gases, because 
(i) the boundary and gravity effects are strong in the rapid granular flows, (ii) the hydrody- 
namic interaction between particles in gas-solid mixtures are complicated, and (iii) a dense 
cluster appears in the vibrating experiment. Therefore, it is difficult to achieve free-cooling 
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gases in experiments. 15 

In this letter, to remove such problems, at first, we propose an experimental accessible 
situation to produce granular gases. Second, we demonstrate that VDFs in both dense and di- 
lute granular gases under vertical vibration can obey exponential-like functions when Coulomb 
friction is important through our simulation based on the distinct element method (DEM). 
This exponential VDF disappears immediately after vibration has stopped, i.e., the grains are 
in a free-cooling process. Third, we introduce a phenomenological Langevin equation eq. (5) to 
describe the motion of particles and explain the mechanism that causes the exponential VDF 
to appeas. The detailed analysis of experimental accessible granular gases will be reported 
elsewhere. 

We use a three-dimensional DEM for monodispersed spherical particles. 26 Instead of using 
the Hertzian contact force, we adopt the linear spring model to represent the repulsion of 
contacted spheres. We also include the rotation of spheres and Coulomb slip for tangential 
contact. The equations of motion of the sphere with the diameter d, the mass m and the 
momentum of inertia / = md 2 / 10 at contact are given by 

mxj = ^2 {Fn^ij + Ft J tij} - mg, 

<ij> 

I* = ^fKxt y }, (1) 

<ij> 

where n^, ty, Yl<ij> an d g represent the normal and the tangential unit vectors at the 
contact point of the j and i particles, the summation over j particles which are in contact 
with i particle, and the gravitational acceleration, respectively. The contact forces F n lJ and 
Ft 3 are respectively given by 

F n %j = -k n (d - |x, - Xj|) - mriniiij • (xj - x,), 

f F t v (if \F t ^\<fi\F n ^\) 
1 fiF n lJ (otherwise) 



F t 13 



(2) 



where 



Ft 3 = -h(d - \xi - Xj|) - mrjttij ■ (x, - x,- + |(ojj + u>j)). (3) 

Here, the friction coefficient /j, describing Coulomb slip is assumed to be fi = 0.8. In addition, 
we introduce four parameters k n , kt,i] n and rjt which are assumed to be k n = 5.0 x W 4 mg/d, 
k t = 2.5 x 10 4 mg/d, rj n = 4.47 x lOy 7 g/d and r\ t = rj n /2 in our simulation. This setup 
corresponds to the normal restitution coefficient e = 0.638. We adopt the second-order Adams- 
Bashforth method for the time- integration with the time interval At = 1.0 x lO^^^d/g. 

We focus on the following situation. Particles are confined in a quasi-two-dimensional 
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Fig. 1. Schematic side view (left) and top view (right) of our setup. The random scatters are fixed 
spherical particles on the top board. The diameter of fixed particles on the top plate (white ones) 
is uniformly random between 0.6c? and 0.8<i, and their centers are located O.lbd above the top 
board. To avoid the stacking of particles in the corners, we introduce four fixed particles (black 
ones in the figure) at the corners. 



container in which the height is 1.8d and the horizontal plane is a square (Fig. 1). In this 
system, particles cannot from a multilayer configuration in the vertical direction, but the 
particles have vertical velocities in their motion. The mobile particles are randomly scattered 
by fixed particles on the top of the container if the container is vertically vibrated. The 
area fractions of fixed scatters are 0.21, 0.23 and 0.25 for 1,000, 3,000 and 10,000 mobile 
particles, respectively. The vibration is driven by a sinusoidal force whose amplitude and 
frequency are given by A = 1.2d and / = 0.5y^gjd. Thus, the acceleration amplitude becomes 
r = A(2irf ) 2 /g = 11.5. If the vibration is stopped, the particles move with the rotation on 
the bottom plate, then the system can simulate a two-dimensional free-cooling granular gas. 
We emphasize that the rotational particles on a flat plane can be regarded as an approximate 
assembly of gas particles in a two-dimensional system, because the energy loss in the rolling 
friction is very small in many cases. 

The simulation starts from a stationary state of particles on the bottom plate. The particles 
gain the kinetic energy from the external oscillation and the system reaches a 'steady state' 
in the balance between collisional dissipations and the gain of the energy from the external 
force. Typically, the system reaches the statistical 'steady state' after 25 oscillation cycles in 
which hydrodynamic quantities are independent of how many cycles proceed from the initial 
condition but depend only on the phase of oscillation. 

In the 'steady state', we confirm that the density is uniform and there is no long-range 
correlation in contrast with granular gases in free-cooling states. 27 Thus, the system does not 
have any systematic flows nor any definite clusters. We also check that the effects of side 
boundaries can be neglected. This is because particles always hit the top wall or the bottom 
wall during the vibration. After we stop the vertical oscillation, the correlation grows with 
time as the free-cooling process proceeds. 

The most significant quantity for charactering this gas system is VDF. In the 'steady 
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Fig. 2. Scaled VDF /(c) in 'steady states' under the vertical vibration. Here, <J/ represents the area 
fraction. For = 7.2% we use 10,000 particles, and use 3,000 particles for other cases. 



Table I. Flatnesses(FN) under several conditions. See text for details. 





Number of 


Fractions(%) 


r IN 




particles 


Area 


Volume 




steady state 


1.0 x 10 3 


6.5 


2.4 


8.00 






10 


3.7 


7.63 






20 


7.4 


6.90 






30 


11 


6.26 






40 


15 


5.67 






50 


19 


5.36 




3.0 x 10 3 


6.0 


2.2 


9.57 






10 


3.7 


8.13 






20 


7.4 


7.03 






30 


11 


6.34 






40 


14 


5.64 






50 


18 


5.27 




1.0 x 10 4 


7.2 


2.6 


6.85 


free-cooling 


1.0 x 10 4 


7.2 


2.6 


4.20 


no tangential 


1.0 x 10 3 


6.5 


2.4 


3.26 




3.0 x 10 3 


6.0 


2.2 


3.43 


no friction 


1.0 x 10 3 


6.5 


2.4 


3.27 




3.0 x 10 3 


6.0 


2.2 


3.43 


undulation 


1.3 x 10 3 






4.10 
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state', the vertical component of VDF has double peaks where each peak corresponds to a 
lifting process or a falling process, while the horizontal VDF has a single peak. For later 
discussion, we are only interested in the horizontal VDFs for analysis. We plot the scaled 
horizontal VDF /(c) in Fig. 2 using 

f(v,t)=nvo(t)- 2 f(v/v ) (4) 

with the density n, the average speed vq = ^2T/m with the granular temperature T, and 
f dcf(c) = f dec 2 /(c) = 1, where /(c) is averaged over the cycles. As in Fig. 2, the scaled 
VDF can be approximated by an exponential function. In fact, the flatness defined by < c 4 > 
/ < c 2 > 2 =< c 4 > with < c n >= J dcc n f(c) is not far from 6. It should be noted that 
the flatness with the Gaussian VDF is 3 and that with the exponential VDF is 6. In our 
simulation, the flatnesses are summarized in Table I for 1,000- and 3,000-particle simulations 
as functions of area fractions projected into the horizontal plane. In addition, the dependence 
of the flatness on the density is relatively weaker in our situation than that by Murayama 
and Sano. 13 A large system with 10,000 particles with an area fraction of 7.2% has a smaller 
flatness of 6.85. If the external oscillation stops, the flatness decreases quickly and is saturated 
at 4.20 for 10,000 particles. As can be seen in Fig. 3, VDF in the cooling process is nearly 
Gaussian for low-energy particles but has an exponential tail for high-energy particles as in 
the usual gas systems. 5-7 

Let us consider the origin of the exponential-like VDF. It is easy to verify that the ex- 
ponential VDF cannot be reproduced without the existence of Coulomb friction in eq. (2). 
Namely, if we eliminate the tangential component of the contact force or the slip rule for 
\F t ^\ < n\F n ^\, VDF becomes a Gaussian-like distribution (Fig. 3). In fact, if all the com- 
ponents of tangential contact force are omitted, the flatness becomes 3.43 for 3,000 particles 
with the area fraction of 6.0%. While if we only neglect the effect of Coulomb friction in 
eq. (2), the flatness becomes 3.43 in the same situation (Table I). In our system, particles 
experience an strong shear force when they hit the fixed scatters on the top wall, and their 
direction of motion changes drastically. Thus, the tangential slips between particles and the 
fixed scatters are the dominant dissipative processes in the 'steady state'. Since we specify 
the origin of exponential-like VDF, we can understand the weak dependence of VDF on the 
density. Namely, particles directly hit the fixed scatters for dilute case, while dense particles 
collide with each other and rotate without slips besides the collisions with the walls. 

Therefore, the essence to produce the large flatness in VDF is apparently Coulomb fric- 
tion. Thus, we propose the following Langevin equation to describe the horizontal motion of 
particles. 

<^ u u / f\ 

= - 7 _-V$ + T7 5 
dt u 

Here, u, u = |u[ and are respectively the velocity, the speed and the potential exerted 



5/9 




Fig. 3. Scaled VDF for several situations. In the inset, 'no tangential' and 'no friction' correspond to 
the results of simulations without the tangential force and the friction force, respectively, where 
the systems are in the 'steady state' and include 3,000 particles with the area fraction 6.0%. Here, 
'cooling state' and 'undulation' indicate the results of the simulation of the cooling process of 
10,000 particles and the undulation, respectively. 



among particles. The friction coefficient 7 may be proportional to fig. The a component of 
the random force rj satisfies the fluctuation-dissipation relation. 

< r) a (t) >= 0, < Va{t)Vf3(t') >= 2D6 a p6(t - t') (6) 



Here, D = 7-y/T '/3m is the diffusion coefficient in the velocity space. The Langevin equation 
eq. (5) with eq. (6) can be converted into the Fokker-Planck or Kramers equation for the 
probability distribution function P(x, u, t). 

at m au 

For spatial homogeneous 'steady states', the equation for P — ► P st is reduced to (u/u)P s t + 
^T/3m(d/du)P st = for u / 0. Its solution can be obtained easily 



P st (u)=2^exp[-^u]. (8) 

Thus, the Langevin equation with Colomb friction law obeys the exponential VDF. 

To check the validity of the new Langevin equation for the motion of particles, we evaluate 
both the diffusion coefficient D and the friction coefficient 7. At first, we have confirmed that 
the motion of particles in the horizontal plane is diffusive. Then, the diffusion coefficient is 
evaluated from the relation < (r(t) — r(io)) 2 >= 4Z)(t — to), where we choose to as 25 cycles 
of the oscillation and simulate the motion of particles until 75 cycles. For each parameter 
of the oscillation we determine T from the second moment of VDF. Thus, we obtain the 
relation between D and T for simulations of 1000 and 3000 particles as in Fig. 4. The best 
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fitted relation is D = 0.0735T 452 - 0.00 78, where D = D /(d 1 / 2 g^ 2 ) and f = T/(mgd), 
but this can be replaced by D = 0.0744 x/t — 0.0054. This result seems to be independent of 
the system size. We also note that the collective vertical motion suppresses the diffusion for 
larger T, and no diffusion takes place because of the insufficient kinetic energy for smaller T. 
From the relation between D and T, we evaluate 7 as 7 = 0.129g. Note that 7 is independent 
of T, because the dominant dissipative process involves collisions between particles and the 
horizontal walls caused by the collective motion of particles in the vertical direction. Thus, 
we believe that the Langevin equation eq. (5) can be used to describe the motion of particles. 

Figure 3 also contains the VDF obtained from the simulation of granular particles confined 
in a dense cluster which exhibits undulations. In the simulation, we use k n = 2.0 x W 4 mg/d, 
h = k n /2, rj n = r/ t = 3.8 x l0^/g~/d and [i = 0.4 with polydispersed particles whose diameters 
are uniformly distributed between 0.8<i and d. The number of particles is 1,300 and confined 
in a thin box whose depth is 4.5<i and width is 60d. The initial accumulated height of particles 
is approximately 6 layers. In this situation, the relative motion of particles is suppressed and 
particles are confined in a lattice-like cage. In this case, VDF is deviated from the exponen- 
tial function where the flatness is approximately 4.10(Table I), but the high-energy particles 
apparently obey an exponential law as in the free-cooling case(Fig. 3). It is easy to imagine 
that particles with a high energy can slip on the cage and particles with a low energy oscillate 
within cages. This picture may be valid for the gas system in the free-cooling process. Thus, 
VDF can be approximated by the combination of the Gaussian part for low-energy particles 
and the exponential part for high-energy particles. The VDFs for the undulation and the free- 
cooling process are similar to those in experiments under vertical vibration in which a dense 
cluster exists and each VDF can be approximated by a single stretched exponential func- 
tion. 11 ' 12,14-18 We consider, however, that the stretched exponential VDF can be understood 
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by the combination of the Gaussian part and the exponential part as stated here. 

We also recall that VDF obeys a Gaussian-like function for dense flows on an inclined 
slope. 28 This is because particles can move without slips because the lattice-like cage is very 
weak. We also compare our results with those of Murayama and Sano. 13 Because of the 
nonexistence of random scatters in their simulation, the slip processes can occur as results of 
collisions among particles. Thus it is reasonable to get the transition from the Gaussian to 
the exponential as the density increases. On the other hand, in our case, the slips take place 
mainly in collisions between particles and fixed random scatters. Therefore the tendency to 
obey the exponential is emphasized for dilute gases. 

This mechanism to obtain the exponential distribution function may be used for other 
situations, because Coulomb friction plays roles in sliding frictions. 29 To look for possibilities 
of the universal feature of our scenario will be our future task. 

In conclusion, we propose an experimental accessible system for granular gases. The VDF 
in a 'steady state' obeys an exponential-like function but changes to Gaussian-like distribution 
function when free-cooling starts. This exponential VDF is caused by Coulomb friction force. 
Thus, we propose the Langevin equation with Coulomb friction to reproduce the results of 
our simulation. 
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